#include <iostream>//must be removed
#include <fstream>
#include <string>
#include <cstdlib>
#include <math.h> //must be removed

#include "excepts.h"
#include "geometry.h"


#include "mesh.h"
#include "globals.h"
#include "neuread.h"
#include "prop.h"
#include "exports.h"

using namespace std;
using namespace ABAN;
//----------------------------------------------------------------------------
/*void showFractCoord(CellIterator X)
{
  cout<<"Fractional Coordination for cell no: "<<X->no<<endl;
  for (int i=0;i<3;i++)
  {
    cout<<X->fractCoord[i].numerator<<"/"<<X->fractCoord[i].denominator;
    cout<<endl;
  }
}*/
//----------------------------------------------------------------------------
/*void showNeighbors(CellIterator X)
{
  cout<<"Neighborhood for cell no: "<<X->no<<endl;
  for (SeqType i=dirXP;i<noOfDirections;i++)
    {
      cout<<"Neighbor "<<i<<" "<<X->neighbors[i]->size()<<" ";
      if (X->neighbors[i]->size())
       for (CellNeighborhoodSetIterator Y=X->neighbors[i]->begin();Y!=X->neighbors[i]->end();Y++)
        cout<<(*Y)->no<<",";
      else
        cout<<"NONE ";
      cout<<endl;
    }

  cout<<endl;

}*/
//----------------------------------------------------------------------------
int main(int argc, char *argv[])
{

  nodes = new VectorContainer();
  cells = new CellContainer();

  cout<<"Reading mesh files... "<<flush;
  ReadNEU::run(nodes,cells,"./mesh/10x10x10.neu"
                          ,"./mesh/10x10x10.abanboundary");
  cout<<"done."<<endl;

  cout<<"Processing neighborhoods... "<<char(13)<<flush;
  cells->checkForMatchingNeighbors();
  cout<<char(13)<<"Processing neighborhoods... "<<flush;
  cout<<"done."<<endl;

  delete nodes;

   //Printing Neighborhood
   CellIterator a,b,c;
   a = b = c = cells->begin();
   for (SeqType i = 0; i<35;i++) a++;
   c = b = a;
   b++;

   cells->breakCell(a,dirXP);
//   cells->breakCell(a,dirYP);
//   cells->breakCell(b,dirYP);
//   cells->breakCell(b,dirYP);
//   cells->breakCell(b,dirYP);
//   cells->breakCell(a,dirYP);
//
//   cout<<"-------------------------------"<<endl;
// //   showFractCoord(b);
//
//   showNeighbors(a);
//   c = cells->end();
//   c--;c--;
//   cells->coarseCell(b,c,dirYP);
//   c = cells->end();
//   c--;c--;
//   cells->coarseCell(b,c,dirYP);
//   showNeighbors(b);
//
//   delete nodes;
//   delete cells;


//  CellIterator a;
//  a=cells->begin();
  //for (SeqType i = 0; i<435;i++) a++;
//  cells->breakCell(a,dirXP);
//  cells->breakCell(a,dirXP);
//  cells->breakCell(a,dirXP);
/*  cells->breakCell(a,dirXP);
  cells->breakCell(a,dirYP);
  cells->breakCell(a,dirXP);
  cells->breakCell(a,dirZP);
  cells->breakCell(a,dirZP);
  cells->breakCell(a,dirZP);*/

  for (CellIterator toBeTested = cells->begin(); toBeTested!=cells->end(); toBeTested++)
  {
    cells->vels[dirXP]->push_back( toBeTested->center[dirXP]*1.0+
                         toBeTested->center[dirYP]*2.0+
                         toBeTested->center[dirZP]*0.0+1.0);
//    cells->vels[dirXP]->push_back( toBeTested->center[dirXP]*0.0+
//                         toBeTested->center[dirYP]*toBeTested->center[dirYP]*(-6.0)+
//                         toBeTested->center[dirZP]*0.0+1.5);
    cells->vels[dirYP]->push_back( toBeTested->center[dirXP]*0.0+
                         toBeTested->center[dirYP]*0.0+
                         toBeTested->center[dirZP]*0.0);
    cells->vels[dirZP]->push_back( toBeTested->center[dirXP]*0.0+
                         toBeTested->center[dirYP]*0.0+
                         toBeTested->center[dirZP]*0.0);
    cells->p->push_back( toBeTested->center[dirXP]*0.0+
                         toBeTested->center[dirYP]*0.0+
                         toBeTested->center[dirZP]*0.0);
    cells->pGradients[dirXP]->push_back(0.0);
    cells->pGradients[dirYP]->push_back(0.0);
    cells->pGradients[dirZP]->push_back(0.0);
  }


//
//    for (SeqType i = 0; i<cells->size();i++)
//        cout<<"cell no "<<i<<"    "<<(*cells->vels[dirXP])[i]<<"    "<<(*cells->vels[dirYP])[i]<<"    "<<(*cells->vels[dirZP])[i]<<endl;

    VTKExport::run(cells,"./out/startcond.vtk");
//    EulerBackwardPropellant* ebp;
//    ebp=new EulerBackwardPropellant(1e-2,1.0,1e-2);
//    ebp->calc(cells);

//    for (SeqType i = 0; i<cells->size();i++)
//        cout<<"cell no "<<i<<"    "<<(*cells->vels[dirXP])[i]<<"    "<<(*cells->vels[dirYP])[i]<<"    "<<(*cells->vels[dirZP])[i]<<endl;


    return EXIT_SUCCESS;
}
